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ARO  Report 

Tools  for  Condensed  Phase  Computational  Chemistry 

An  Army/University  of  Delaware  Workshop  held  at 
Clayton  Hall,  Newark,  Delaware 
27-30  May  2003 

Organizers: 

Prof.  Krzysztof  Szalewicz,  Department  of  Physics,  University  of  Delaware 
Dr.  Robert  W.  Shaw,  Chemical  Sciences,  Army  Research  Office 


Abstract: 

This  workshop  brought  international  experts  on  condensed  phase  computational  chemistry  codes 
together  with  Army  scientists  working  on  energetic  materials  to  help  ensure  that  Army  scientists 
have  available  the  best  current  computational  tools  (codes).  The  workshop  was  planned  by  U 
Delaware,  ARO,  and  ARL-WMRD. 

This  workshop  report  is  a  guide  for  the  Army  research  program  in  energetic  material  design. 
New  research  collaborations  between  Army  scientists  and  experts  at  the  workshop  are  likely  and 
ongoing  collaborations  Were  improved. 
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Description  of  the  Workshop:  Tools  for  Condensed  Phase  Computational  Chemistry 

The  workshop  developed  recommendations  for  those  computational  chemistry  tools  best  suited 
to  attack  research  problems  in  properties  of  energetic  molecules  -  individually  and  in  condensed 
phases.  A  small,  international  group  of  experts  presented  their  work  and  discussed  current  tools 
of  computational  chemistry  in  the  context  of  the  Army  program  of  energetic  material  design. 
These  experts  (see  Participants  List),  chosen  for  their  specific  skills,  came  from  Europe  and  the 
US. 

The  workshop  emphasized  problems  that  can  be  treated  now.  To  facilitate  discussion,  we 
described  the  tools  as  falling  into  three  groups: 

1 .  Readily  available,  widely  used,  relatively  fnendly  packages,  e.g.,  MolPac,  Gaussian, 

Cerius  suite 

2.  Available,  but  directed  to  the  specialist,  require  special  training,  e.g.,  DMAREL 

3.  Written  and  used  within  a  research  group,  accessible  only  to  a  few. 

Questions  raised  included:  What  tools  are  best  suited  to  current  Army  computational  problems? 
What  investment  (time,  money)  is  required  to  use  these  tools?  How  can  the  tools  output  be 
validated  (what  experiments)? 

The  workshop  began  with  Army  scientists  (Betsy  Rice,  Cary  Chabalowski,  William  Mattson) 
describing  their  current  goals  and  problems.  These  included  predictions  of  crystal  structures  and 
other  physical  properties  of  energetic  molecules.  Visiting  experts  came  from  academia 
(University  College-London,  Cambridge  University,  Imperial  College,  Rice  University,  UCLA, 
Mississippi  State  University,  and  Universities  of  Illinois,  Delaware,  and  Maryland),  Naval 
Research  Lab,  Pacific  Northwest  Lab,  and  NIST.  These  visiting  experts  described  capabilities  of 
current  computational  chemistry  codes. 

All  these  presentations  included  time  for  questions  and  discussion.  Following  the  presentations, 
the  participants  divided  into  groups  led  by  Army  scientists  to  discuss  which  packages  would  best 
suit  current  Army  projects  and  what  research  should  be  done  to  develop  these  packages  to  enable 
future  Army  projects  in  design  of  energetic  materials.  At  the  end,  these  discussions  were 
summarized  and  discussed  by  everyone.  This  report  includes  those  summaries. 

The  workshop  scientific  sessions  and  meals  were  held  at  the  Clayton  Hall  Conference  Center. 
Participants  all  lodged  in  the  same  hotel  (Embassy  Suites)  to  enhance  opportimities  for 
discussion. 

Previous  Army  Statements  of  Interest 

The  discussions  were  based,  in  part,  on  questions  raised  at  the  ARO/ARL/ARDEC  200G  EM 
Planning  Meeting.  Some  specific  questions  follow. 

Computational  Chemistry  and  Properties  of  Energetic  Molecules 
Question:  What  theoretical  tools  are  most  important? 

Answer:  We  need  reliable  force  fields  (many  body  forces)  to  model  chemical  reactions  in  the 
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condensed  phase. 

Computation  requirements  for  the  eondensed  phase  are  large  beeause  we  must  inelude  many 
moleeules  (many  body  problem).  Often  these  condensed  phase  ealeulations  rely  on  periodieity  in 
the  structure  to  reduee  computation  time.  But  chemical  reactions  break  this  order. 

We  should  move  our  methods  to  determine  low-energy  crystal  structures  of  energetie  materials 
from  elassieal  molecular  packing  to  quantum  mechanical  molecular  packing.  This  proj  ect  may 
be  appropriate  for  a  graduate  student  thesis. 

Sereening  Tools:  We  want  to  use  theory  to  screen  notional  materials  (moleeules  that  may  not 
have  been  synthesized  but  whieh,  based  on  their  expected  structures,  appear  to  have  good 
properties).  The  first  properties  we  want  to  predict  are:  burning  rate,  heats  of  formation, 
densities. 

Density:  Moleeular  packing  calculations  are  the  best  method  for  predicting  densities  of  new 
solid  materials.  One  must  generate  a  variety  of  crystals  corresponding  to  the  various  erystalline 
spaee  groups  from  the  coordinates  of  a  single  molecule.  Each  crystal  must  be  optimized  (the 
erystallographic  parameters  are  minimized  with  respect  to  energy,  while  the  molecular 
parameters  themselves  remain  fixed)  and  all  energies  are  compared.  The  current  criterion  for 
seleetion  of  the  most  probable  structure  is  lattice  energy;  i.e.,  the  structure  that  has  the  lowest 
lattiee  energy  is  assumed  to  be  the  “correct”  crystal  structure.  This  criterion  has  proven  to  be 
insufficient  and  additional  criteria,  perhaps  the  evaluation  of  AG  for  each  crystal  structure,  should 
be  included  in  the  selection  procedure. 

Additional  improvements  to  such  calculations  would  be  improved  intermolecular  interaction 
potentials,  and  relaxation  of  the  rigid  molecule  assumption  to  partial  deformation  of  the 
molecules  (inclusion  of  torsional  rotations  for  example).  Although  molecules  like  RDX  may  be 
adequately  modeled  as  rigid  structures,  allowance  for  molecular  deformation  will  likely  be 
especially  important  for  floppy  molecules  like  PETN.  So  we  should  consider  a  stepwise 
approach:  optimize  AG  (requiring  calculation  of  the  energy  second  derivatives  of  the  erystals, 
very  time  consuming)  and  allow  some  intramolecular  motion. 

All  calculations  discussed  to  date  are  classical;  hence  they  require  an  intermolecular  force  field. 
And,  as  mentioned  above,  do  not  include  intra-molecular  dynamics.  The  energetie  moleeule 
intramoleeular  potential  developed  by  Sorescu,  Rice,  and  Thompson  (1)  has  suceessfully 
matehed  experimental  crystallographic  structures  and  lattice  energies  for  a  variety  of  explosives. 
However,  the  foree  field  was  tested  by  generating  the  crystal  structure  of  RDX  assuming  the 
experimentally  determined  space  group,  and  doing  an  energy  minimization.  These  results  give 
erystal  structures  in  good  agreement  with  experiment.  Thus,  predietions  of  densities  are  good.  It 
may  be,  however,  that  the  force  field,  in  the  absence  of  space  group  information,  would  prediet 
other  erystal  structures  with  more  favorable  AG.  To  truly  test  the  foree  field,  one  must  paek 
RDX,  for  example,  in  every  possible  fashion  according  to  the  230  erystalline  spaee  groups  and 
find  the  best  AG.  To  have  eonfidence  in  the  force  field,  it  must  generate  the  eorreet  erystalline 
spaee  group  without  relianee  on  experiment  at  all. 


1  J.  Phys  Chem  A,  103,  989-998,  1999 
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Summary  by  Focus  Group  1 

(reported  by  Cary  Chabalowski,  Ed  Byrd,  William  Mattson) 

The  discussion  session  was  divided  into  two  general  topics; 

What  do  we  need  to  do  to  predict  sufficiently  accurate  crystal  densities  and  other  static 
properties?  (At  least  better  than  what  is  currently  available). 

What  do  we  need  to  predict  the  chemistry  of  detonation  in  the  solid  state?  Along  with  this,  we 
asked  that  both  short-term  and  long-term  goals  be  put  forth  in  terms  of  modeling  developments. 

Comments: 

When  we  refer  to  optimizing  crystal  structmes,  we  mean  the  ability  to  optimize  both  cell 
parameters  and  atomic  positions.  PBC  stands  for  periodic  boundary  conditions. 

TOPIC  #1:  Crystal  Densities 

This  question  seemed  to  evolve  into  discussions  on  methods  for  predicting  the  complete  van  der 
Waals  forces.  There  was  also  discussion  on  sampling  configuration  space.  Some  of  the  recent  work 
by  Mark  Tuckerman  may  be  of  interest  here  (conformational  sampling  algorithms). 

l.A  Force-Field  Approach 

The  currently  used  methods  of  predicting  crystal  structure  are  based  on  the  use  of  (mostly 
empirical)  force  fields.  The  discussion  concentrated  on  possible  applications  of  non-empirical 
methods  to  develop  such  force  fields. 

Asymptotic  force  fields: 

Not  many  new  ideas  were  proposed  on  the  topic  of  predicting  accurate  electrostatic  intermolecular 
interactions,  since  a  general  sense  prevailed  that  currently  existing,  or  soon-to-be  developed.  Density 
Functional  Theory  (DFT)  functionals  will  be  sufficiently  accurate  to  give  good  asymptotic 
electrostatic,  induction  (polarization),  and  dispersion  interaction  energies.  The  latter  may  sound  as 
going  against  the  common  knowledge,  but  is  true  if  time-dependent  DFT  (TD-DFT)  is  used  to 
calculate  the  dynamic  polarizabilities,  see  below.  However,  some  discussion  did  occur. 

Force-fields  including  overlap  and  exchange  effects: 

In  order  to  calculate  the  force  fields  in  the  region  of  the  van  der  Waals  minimum,  critical  for 
determining  the  crystal  structure,  one  has  to  include  the  overlap  (charge  penetration)  and  exchange 
effects.  The  overlap  effects  significantly  change  the  values  of  the  electrostatic,  induction,  and 
dispersion  forces  in  this  region  and  the  exchange  energies,  providing  the  repulsive  forces,  are 
responsible  for  the  very  existence  of  van  der  Waals  minima. 
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(LONG-TERM  GOAL  for  supermolecular  DFT  methods) 

Due  to  its  efficiency,  the  DFT  method  would  potentially  enable  calculations  of  interaction  force 
fields  for  molecules  typical  for  energetic  materials.  Unfortunately,  current  exchange-correlation 
DFT  functionals  have  shown  mixed  results  in  predicting  accurate  electrostatic  and  induction 
interactions  with  charge-penetration  effects  (as  discussed  in  Szalewicz's  talk),  including  hydrogen 
bonding.  It  is  believed  to  be  due  to  the  fact  that  asymptotic  corrections  (AC),  fixing  up  the  electron 
density  at  large  distances  from  the  center  of  a  molecule,  have  not  been  applied  (see  below).  Actually, 
it  is  not  even  clear  how  to  apply  these  corrections  in  supermolecular  approach.  In  any  case,  a  more 
serious  failure  of  DFT  is  in  recovering  the  dispersion  energy  (see  below).  However,  the  asymptotic 
electrostatic  and  induction  energies,  depending  only  on  molecules'  multipole  moments  and  static 
polarizabilities,  are  predicted  well  by  the  modem  DFT  even  without  AC.  Thus,  for  small  and  very 
polar  molecules,  whose  interactions  are  dominated  by  the  asymptotic  electrostatics  and  induction, 
supermolecular  DFT  approach  may  give  reasonable  predictions.  Such  predictions  will  be  further 
improved  by  ongoing  research  into  the  development  of  better  DFT  functionals,  independent  of  any 
specific  work  on  developing  better  solid  state  modeling  methods.  Using  supermoleculm*  DFT  to  deal 
with  general-case  van  der  Waals  forces  is  totally  an  open  issue  and  an  unsolved  problem.  Several 
groups  investigate  this  problem  but  it  is  impossible  to  say  when  and  if  a  solution  will  be  found. 

(SHORT-TERM  GOALS  for  building  non-empirical  force  fields) 

One  possible  method  for  obtaining  accurate  electrostatic,  induction,  and  exchange  interaction 
energies  valid  at  all  intermolecular  separations  is  through  the  use  of  a  newly  developed  version  of 
Jeziorski-Szalewicz’s  Symmetry  Adapted  Perturbation  Theory  (SAPT),  which  is  based  upon  the  use 
of  the  Kohn-Sham  (KS)  orbitals  from  noninteracting  monomers  in  the  SAPT  expansion,  as  proposed 
by  Williams  and  Chabalowski.  This  approach,  referred  to  as  SAPT(KS),  when  used  with 
asymptotically  corrected  functionals,  gave  electrostatic,  induction,  and  exchange  interaction  energies 
for  dimers  that  compared  to  within  a  few  percent  of  benchmark  traditional  SAPT  calculations,  but  at 
a  much  reduced  computational  cost.  The  size  of  the  monomers  that  can  be  treated  is  yet  to  be 
determined,  but  Szalewicz  thought  that  a  dimer  of  RDX  could  be  studied  in  near  future.  The 
remaining  component  of  the  interaction  energy,  the  dispersion  force,  can  also  be  accurately  computed 
using  SAPT  based  on  TD-DFT  dynamic  polarizabilities.  These  dimer  interaction  energies  can  be 
used  in  a  fitting  procedure  to  create  non-empirical  force  fields  to  be  used  in  modeling  of  molecular 
crystals.  These  energies  can  also  be  used  to  improve  the  individual  potential  energy  terms  in 
anal3dical  empirical  force  fields,  as  discussed  in  the  next  subsection. 

Empirical  Force  Fields: 

(SHORT  TERM  GOAL) 

If  reliable  fitting  data  is  available  through  either  experiment  or  theory  [SAPT(KS)  was  specifically 
mentioned],  it  was  suggested  that  classical  force  field  packages  (e.g.  DMAREL  from  Price),  could 
and  would  be  improved  by  refitting  using  fitting  codes  such  as  GULP.  There  was  a  general  sense 
that  these  force  field  approaches  were  capable,  in  general,  of  reproducing  the  electrostatic 


5 


interactions  in  crystals  to  an  acceptable  accuracy  provided  that  accurate  multipole  moments  are 
known.  Some  discussion  ensued  about  the  need  for  including  an  effective  anisotropic  partial  charge 
distribution,  or  an  equivalent  multipole  expansion  of  the  charge,  around  atoms  or  groups,  allowing 
for  interactions  such  as  those  involving  an  electron  lone-pair.  The  need  to  include  other  refinements 
was  also  discussed,  such  as  polarizable  charge  distributions  or  multipole  expansions  of  the  charge. 
However,  these  refinements  are  probably  less  important  than  the  inclusion  of  a  better  representation 
of  the  remaining  components  of  the  interaction  energy  like  dispersion  or  exchange. 

In  this  spirit,  Murdachaew  suggested  a  specific  approach  to  developing  a  force  field  including 
dispersion  and  other  interaction  terms  based  on  Szalewicz’s  SAPT.  The  ab  initio  SAPT  energies  can 
be  used  to  produce  first-principles  force  fields.  These  energies  are  fed  directly  into  a  fitting  code 
developed  by  Robert  Bukowski  (in  Szalewicz’s  group).  This  code  uses  asymptotic  van  der  Waals 
constants  computed  ab  initio  with  the  POLCOR  code  developed  in  Paul  Wormer’s  group  in 
Nijmegen).  The  resulting  force  field  produces  force  fields  with  accurate  asymptotic  expansions 
expansions  that  include  the  terms  containing  the  Ce  through  Cn  coefficients.  These  coefficients  or 
the  monomer  properties  imderlying  them  can  be  used  to  improve  empirical  force  fields.  The 
monomer  properties  can  also  be  computed. 

Krzysztof  Szalewicz  comments:  What  Garold  Murdacheaw  proposed  makes  a  lot  of  sense.  In  fact,  I 
had  a  long  discussion  with  Sally  Price  after  the  conference  and  I  suggested  to  her  something  similar 
(not  knowing  about  Garold's  suggestion).  They  should  be  able  to  investigate  quite  large  molecules 
using  a  much  better  asymptotics  than  they  do  now,  basically  following  our  approach  to  asyniptotics. 
Notice  that  in  this  way  the  asymptotic  dispersion  interaction,  which  is  included  only  in  a  very  crude 
way  in  current  empirical  force  fields,  will  be  very  well  accounted  for.  Furthermore,  they  could 
improve  the  exchange  by  utilizing  the  approximate  proportionality  of  exchange  to  overlap  effects. 
We  are  now  using  this  approach  for  flexible-monomer  potentials,  but  Sally  Price  with  Richard 
Wheatley  actually  pioneered  it. 

(LONG-TERM  GOAL) 

It  was  suggested  that  a  large  database  (primarily  calculated?)  of  intermolecular  interaction  energies, 
as  well  as  other  molecular  properties,  be  established.  Such  data  could  be  very  useful  in  refining 
analytical  force  fields.  The  current  database  at  NIST  is  an  obvious  effort  in  this  direction,  but  for 
individual  niolecules,  and  should  be  expanded  to  include  inter-molecular  interaction  energies  versus 
structures.  This  data  should  be  assessed  and  documented  for  its  accuracy  whenever  possible,  and 
made  available  to  the  research  community  in  an  easily  searchable  form. 

Dispersion  Forces 

(SHORT-TERM  GOAL) 

One  more  alternative  to  developing  force  fields  for  molecular  crystals  is  to  use  the  Hartree-Fock  SCF 
method  which  gives  a  good  account  of  electrostatic,  induction,  and  exchange  energies  and  at  the 
same  time  is  about  as  fast  as  the  DFT  method.  SCF  does  not  include,  however,  any  dispersion  terms. 
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The  latter  energy  is  a  pure  correlation  effect  and  therefore  its  calculation  will  always  be  much  more 
time  consuming.  One  possibility  is  to  add  the  asymptotic  dispersion  to  the  SCF  interaction  energy, 
but  the  accuracy  of  this  approach  may  not  be  high  enough  for  crystal  modeling.  A  better  accuracy 
will  be  achieved  by  using  the  TD-DFT,  also  called  Coupled  Kohn-Sham  (CKS)  method,  to  calculate 
dynamic  polarizabilites  and  apply  them  in  the  disp(CKS)  method  of  Szalewicz  to  get  dispersion 
interactions.  More  precisely,  this  approach  utilizes  closely  related  quantities.  Frequency  Dependent 
Density  Susceptibilities  (FDDSs)  as  predicted  by  TD-DFT.  The  FDDSs  can  then  be  used  in  a 
generalized  Casimir-Polder  expression  for  dispersion  energy  valid  at  any  intermolecular  separation, 
and  coded  such  that  the  cost  of  computations  scales  as  n^,  where  n  is  the  number  of  molecular 
orbitals  in  a  dimer  calculation.  This  approach,  again  combined  with  asymptotically  corrected 
functionals,  was  capable  of  predicting  dispersion  energies  within  a  few  percent  of  the  benchmark 
SAPT  calculations  for  dimers.  While  this  approach  was  said  to  scale  as  n^,  further  computational 
efficiency  might  be  possible. 

(LONG-TERM  GOAL  for  DFT  dispersion) 

Although  the  supermolecular  DFT  method  is  now  known  to  give  fairly  large  errors  for  the 
electrostatic,  induction,  and  exchange  components,  the  weakest  element  of  this  method  is  the 
dispersion  energy.  This  is  well  known  from  the  fact  that  for  interactions  of  rare  -gas  atoms  which 
are  dominated  by  dispersion,  DFT  gives  simply  nonsensical  results.  Probably  all  the  DFT  exchange- 
correlation  functionals  include  some  portion  of  dispersion  interactions  (for  example,  there  are 
dispersion  interactions  in  the  uniform  electron  gas,  thus,  even  LSDA  contains  some  dispersion),  but 
not  even  the  empirically  fitted  or  hybrid  functionals  include  it  well.  There  have  been  intensive 
efforts  to  develop  functionals  that  could  account  for  dispersion  interactions.  Kohn  in  1998  (and 
perhaps  independently  Lundquist)  proposed  to  use  Time  Dependent  (TD)  DFT  to  calculate  the 
and  Ci2  coefficients  and  include  them  in  the  functional.  This,  however,  constructs  functionals  usable 
only  for  that  specific  molecule  and  not  a  general  functional. 

Szalewicz  comments:  Indeed,  all  the  papers  claiming  to  have  solved  the  DFT  problems  with 
dispersion  interactions  ended  up  in  computing  just  the  asymptotic  C„  dispersion  coefficients. 
However,  this  task  has  been  de  facto  achieved  earlier  by  all  the  DFT  papers  devoted  to  frequency- 
dependent  polarizabilities  since  once  these  are  known,  computation  of  €„  from  the  Casimir-Polder 
asymptotic  formula  is  trivial. 

l.B  Periodic  Boundary  Conditions  Approach 

Virtually  all  computational  predictions  for  molecular  crystals  involving  large  organic  molecules  are 
obtained  with  force-field  methods.  The  PBC  approach  where  optimization  of  the  lattice  can  be 
performed  on  the  potential  energy  surface  computed  “on-the-fly”  would  offer  some  advantages,  in 
particular  where  the  relaxation  of  monomer  structures  is  important.  However,  the  general  consensus 
was  that  no  computationally  feasible  approach  currently  exists  for  predicting  crystal  cell  structures 
using  PBCs  that  give  sufficiently  accurate  predictions  of  intermolecular  dispersion  interactions, 
including  DFT. 
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MP2-PBC  method 

(LONG-TERM  GOAL  for  crystal  relaxations) 

The  MP2  based  codes  of  Scuseria  (the  PBC  code  is  available,  the  MP2-PBC  is  currently  not 
available  in  Gaussian-03)  will  account  for  dispersive  interactions  in  a  PBC  framework  within  the 
limit  of  the  accuracy  of  MP2  theory.  While  MP2-PBC  exists,  it  is  still  a  fledging  code/theoiy,  and 
requires  work  to  fully  develop  and  test.  Unfortunately,  at  this  time,  there  are  no  analytical  gradients 
in  G03  within  the  MP2-PBC  approach,  so  minimum  energy  crystal  structures  cannot  be  readily 
determined.  Scuseria  said  that  he  plans  to  code  these  gradients,  but  this  will  take  time.  Other 
concerns  about  the  use  of  MP2  included  the  inability  to  handle  metalization,  wherein  the  band  gap 
decreases  to  near  zero.  This  is  a  fundamental  weakness  in  the  use  of  MP  theory,  which  requires  the 
use  of  band  gap  energies  (virtual  orbital  energies  in  molecules),  and  it  was  not  clear  how  to  correct 
this  problem.  Another  issue  surrounding  the  use  of  MP2  was  the  computational  efficiency.  In 
principle,  MP2  does  not  scale  as  well  as  DFT  with  respect  to  the  number  of  electrons.  Work  needs  to 
be  done  to  improve  computational  efficiency  to  allow  for  applications  to  larger  unit  or  supercells. 

Emily  Carter  comments:  The  problem  is  that  MP2  conventionally  scales  as  n^  and  all  the  truncation 
tricks  used  by  all  of  us  doing  reduced  scaling  cannot  be  used  here  because  the  whole  point  is  to  retain 
the  long-range  correlation  to  get  the  van  der  Waals  forces  right.  So  it  is  not  clear  at  all  how  to  make 
MP2  scale  in  a  way  that  makes  it  tractable  without  messing  up  the  long-range  interactions  you  want 
to  keep. 

Meanwhile,  Gustavo  Scuseria  states:  Our  MP2-PBC  is  linear  scaling,  it  is  the  prefactor  that  is  huge, 
so  calculations  are  expensive.  Even  the  HF-PBC  is  expensive.  Work  needs  to  be  done  to  improve 
efficiency  and  computer  implementation  of  the  current  approach,  which  is  still  in  its  infant  stages. 

(LONG-TERM  GOAL) 

The  suggestion  was  made  to  develop  DFT  OEPs,  Vxc,  by  optimizing  the  orbitals  within  an  MP2 
calculation.  Apparently,  the  idea  was  to  generate  exchange-correlation  operators  that  retained  the 
capabilities  to  generate  electron  densities  that  reproduce  the  MP2-like  van  der  Waals  forces. 

Emily  Carter  says:  As  far  as  I  understand  the  OEP  strategy,  you  still  have  to  do  the  MP2  calculations 
(an  n^  process  as  we  were  just  discussing)  and  then  one  derives  the  KS  EP  from  that.  And  no  one 
has  figured  out  how  to  go  from  one  of  those  EP’s  to  derive  new  general  functionals,  as  far  as  I  know. 

Corrections  to  force-fields 

(SHORT-TERM  GOAL) 

Professor  Doren  states:  One  method  that  depends  only  on  existing  computational  methodology  is 
a  systematic  correction,  based  on  ab  initio  calculations  to  the  predictions  of  force-field  methods. 
The  calculation  begins  with  a  crystal  supercell,  constructed  and  optimized  with  the  “best 
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available”  force  field  approach.  The  force  field  is  used  to  run  a  dynamics  simulation,  from  which 
snapshots  are  taken  to  provide  an  ensemble  of  configuration  space  samples.  Ab  initio  single- 
point  energies  for  the  sampled  configurations  can  be  calculated  using  either  periodic  MP2  (e.g., 
with  G03),  or  by  calculating  interactions  of  a  unit  cell  with  a  small  number  of  neighbor  cells 
using  supermolecule  MP2  (or  possibly  SAPT)  methods.  Using  free  energy  perturbation,  a 
suitable  ensemble  average  oveftthese  configurations  provides  an  ab  /n/t/o  correction  to  the  free 
energy,  including  dispersion  at  the  level  of  the  ab  initio  calculation.  The  method  has  been 
applied  successfully  to  liquid  state  solvation,  where  50-100  configurations  typically  suffice  for 
accurate  corrections.  Moreover,  in  cases  where  the  ab  initio  calculations  show  that  the  force 
field  model  is  inaccurate,  the  ab  initio  results  can  be  used  to  refit  the  force  field.  Thus,  the 
method  provides  both  a  test  of  the  accuracy  of  force  fields,  and  a  correction  to  their  predictions. 
The  method  can  be  readily  implemented  without  new  technical  developments,  and  given  a 
reasonable  force  field,  the  cost  of  configuration  sampling  is  inherently  lower  than  that  of  using  an 
ab  initio  method  for  sampling. 

Monte-Carlo  methods 

(LONG-TERM  GOAL) 

A  more  rigorous  approach  to  phase  space  sampling  was  suggested  by  Martin  in  the  form  of 
improving  Quantum  Monte  Carlo  (MC)  calculations.  He  expanded  this  idea  by  suggesting  a 
correlated  sampling  combined  with  Diffusion  MC.  He  stated  that  better  correlated  sampling  is 
needed  to  converge  van  der  Waals  forces. 

Embedding  approach 

(LONG-TERM  GOAL) 

The  suggestion  was  made  to  pursue  an  embedding  approach.  A  high-level  QC  method  is  embedded 
into  a  lower-level  treatment,  so  dispersion  is  explicitly  included  within  the  “embedded”  part  of  the 
model  via  QC.  The  remainder  of  the  system  still  faces  the  same  issues  of  accounting,  in  some 
suitable  fashion,  for  the  dispersive  interactions.  Some  questions  were  raised  concerning  whether  or 
not  this  type  of  approach  could  be  used  with  PBCs,  but  the  response  seemed  to  be  that  it  could. 

Emily  Carter  adds:  The  only  reason  to  think  about  doing  this  is  to  see  if  the  dispersion  calculated  via 
MP2  embedded  in  a  periodic  background  treated  with  DFT  (something  we  can  do  today)  gives  a 
more  accurate  description  than  does  MP2  on  isolated  dimers. 

SAPT  with  PBC 

(LONG-TERM  GOAL) 

One  serious  issue  concerning  SAPT  is  whether  or  not  it  can  be  formulated  for  use  with  periodic 
boundary  conditions.  There  had  been  some  behind-the-scenes  discussions  during  the  conference 
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between  Scuseria  and  Szalewicz,  resulting  in  the  notion  that  this  might  be  possible.  After  the 
conference,  Szalewicz,  in  consultation  with  Bogumil  Jeziorski  from  University  of  Warsaw, 
developed  a  strategy  for  such  a  case.  This  strategy  build  upon  the  existing  SAPT  approach  for  trimers 
which  can  be  applied  to  larger  clusters  as  well  without  new  conceptual  difficulties.  The  current 
SAPT  approach  analytically  splits  the  interaction  energy  into  pair  interactions  and  nonadditive  terms. 
This  would  be  difficult  for  larger  clusters,  but  if  the  splitting  is  not  preformed,  the  method  may  be 
practical  and  allows  the  use  of  PBC. 

To  define  this  method  more  precisely,  consider  a  set  of  atoms  or  molecules  in  a  cell.  This  cell  does 
need  to  be  particularly  large,  typical  size  used  in  DFT  PBC  calculations  should  be  sufficient. 
However,  the  cell  should  contain  several  monomers.  If  we  neglect  all  the  terms  in  the  Hamiltonian 
that  describe  interactions  between  the  molecules,  this  defines  our  operator  Ho  and  the  removed  terms 
are  the  operator  V.  We  solve  the  zeroth-order  problem  and  build  the  perturbation  theory  on  top  of  it, 
all  with  PBC.  There  may  be  some  issues,  hopefully  solvable  somehow,  related  to  simultaneously 
imposing  the  permutational  and  translational  symmetry.  Like  in  the  regular  SAPT,  the  zeroth-order 
wave  function  is  not  fully  antisymmetric  under  permutations  of  electrons.  Therefore,  “symmetry 
adaptation”  used  in  the  regular  SAPT  has  to  be  performed  here  as  well.  This  must  be  done  in  a  way 
preserving  translational  symmetry  of  the  crystal. 

SAPT  can  also  be  developed  to  calculate  interactions  of  a  molecule  with  surface.  One  monomer  is 
the  surface,  represented  by  a  few  layers  of  crystal  molecules/atoms  in  a  very  large  supercell.  The 
other  monomer  is  the  single  molecule  in  the  same  supercell.  The  supercell  is  repeated  periodically 
in  all  directions.  Then  one  can  calculate  the  interaction  of  these  two  monomers  using  SAPT. 

An  alternative  way  Szalewicz  can  see  to  extend  it  to  the  bulk  crystal  is  to  assume  a  fairly  large 
supercell  with  one  of  the  molecules  removed  as  one  monomer.  The  other  monomer  would  be  a  single 
molecule  in  the  same  supercell,  placed  in  the  position  of  the  void  in  the  first  monomer.  Again,  PBC 
are  imposed  on  supercell. 

TOPIC  #2:  Chemistry  of  Detonation  in  the  Solid  State 

Predicting  accurate  reaction  pathways  is  still  a  cutting-edge  research  problem  even  for  unimolecular 
reactions,  let  alone  in  crystals.  However  before  discussing  accuracy  of  methods,  a  question  arose 
concerning  the  usefulness  of  minimum  energy  pathways  in  shock  and  detonation  studies  of  energetic 
crystals:  Are  transition  states  for  minimum  energy  paths  useful  in  modeling  shock  and  detonation, 
because  the  input  energy  from  the  shock  wave  might  be  so  great  as  to  raise  the  internal  energy  of  the 
molecules  far  above  the  reaction  threshold.  If  so,  the  atoms  might  not  sense  the  transition  states. 
True?  False?  TBD?  For  the  same  reason,  is  knowledge  of  the  minimum  energy  pathway,  in 
general,  even  useful? 

Emily  Carter  notes:  If  one  does  not  even  characterize  the  pathways  at  zero  K  and  zero  P  first,  I  don’t 
see  how  one  can  have  any  idea  how  to  characterize  pathways  at  arbitrary  T  and  P.  I  think  it  is 
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essential  to  map  out  the  PES  as  much  as  one  can,  as  a  function  of  P,  and  then  one  can  know  how  T 
may  affect  the  pathways  observed  for  a  given  P. 

(SHORT-TERM  GOAL) 

On  a  positive  note,  assuming  that  current  DFT  functionals  could  accurately  predict  chemical 
reactions,  there  are  already  methods  available  to  study  aspects  of  shock  and  detonation.  Solid  state 
codes  have  already  been  used  to  simulate  shock  waves  [we  are  assuming  these  are  pseudopotential 
planewave  codes?].  Shock  waves  have  been  simulated  in  many  ways:  the  continuum  level 
hydrodynamic  approach  (no  input  from  atomistic  level),  from  molecular  dynamics  simulation  using 
analytic  potential  functions,  and  more  recently  by  Carter  and  Ortiz  using  a  multiscale  finite  element 
approach  that  takes  input  from  DFT.  The  current  methods  also  allow  study  of  isotropic  and 
anisotropic  compressions,  both  in  bimolecular  systems  and  in  crystals.  High  temperature  Bom- 
Oppenheimer  Molecular  Dynamics  or  Car-Parinello  MD  simulations  can  also  be  done  (on  a 
reasonable  time  scale  -  one  must  be  careful  here  and  check  the  time  scale  of  the  detonation  wave. 
BOMD  or  CPMD  calculations  can  be  done  for  time  of  the  order  of  ps  but  not  ns)  with  today’s 
methods.  Another  thought  was  to  study  aspects  of  detonation  by  simulating  two  molecules  colliding 
with  large  kinetic  energy. 

Emily  Carter  states:  Take  it  from  someone  who  has  done  simple  bimolecular  collisions,  one  can 
look  at  representative  trajectories  but  you  will  not,  with  ab  initio  MD,  be  able  to  do  enough  statistics 
to  really  sample  phase  space. 

(LONG-TERM  GOAL) 

But  once  again,  the  accuracy  of  DFT  theory,  or  more  appropriately  the  existing  exchange-correlation 
functionals,  was  said  to  be  unreliable  in  predicting  accurate  transition  state  (TS)  structures  and 
energies,  including  the  possibility  of  representing  unrealistic  electronic  spin  states  along  the  reaction 
paths.  What  can  we  do  to  improve  the  functionals  at  the  TS,  and  insure  that  the  densities  represent 
the  correct  spin  state  of  the  system  at  all  points  along  a  reaction  path? 

Emily  Carter  says:  It  is  certainly  true  that  spin  contamination  can  be  a  problem  in  DFT,  imless  one 
does  restricted  DFT,  which  is  not  usually  done.  The  problem  is  that  really  one  needs  a  multireference 
wavefunction  with  all  possible  spin  functions  included  variationally  in  order  to  get  proper  reaction 
paths.  This  brings  us  back  to  ab  initio  methods  however,  since  multireference  DFT  does  not  make 
sense  to  me. 

(UNCLEAR  as  to  state  of  progress) 

Another  problem  with  many  functionals  when  used  in  mapping  reaction  PESs,  is  the  well-known 
electron  “self  interaction  error.”  It  was  suggested  that  methods  already  exist  for  correcting  this  error 
in  some  functionals,  but  the  consensus  seemed  to  be  that  more  work  needs  to  be  done  on  this  topic. 
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Emily  Carter  adds:  Yes  this  is  the  problem  with  radieals.  All  DFT  methods  without  full  exaet 
exehange,  i.e.,  all  of  them,  will  have  this  problem.  More  work  needs  to  be  done  on  this  topie. 

(LONG-TERM  GOAL) 

For  erystals  under  very  high  pressures  or  pressure  waves,  it  was  said  that  a  better  deseription  of  the 
exeited  eleetronic  states  is  needed  within  DFT,  or  at  least  the  exeitation  energies  into  these  states. 
This  problem  with  predieting  aceurate  exeitation  energies  is  also  related,  at  least  in  part,  to  the  well- 
known  problem  of  errors  in  the  eigenvalues  assoeiated  with  the  Kohn-Sham  unfilled  orbitals  (or 
bands)  whieh  tend  to  have  ineorrectly  small  values. 
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Report  of  Focus  Group  2 

(reported  by  Betsy  Rice) 

The  task  of  the  discussion  groups  was  to  determine  what  can  be  done  using  currently  available 
tools  and  to  identify  critical  "deficiencies"  in  the  models/tools  that  deserve  attention  in  order  to 
further  the  methodology.  The  currently  used  methods  of  energy  optimization  with  empirical 
force  fields  cannot  always  identify  the  correct  crystal  structures.  However,  if  the  X-ray  model  is 
used  as  the  packing  probe,  the  smallest  lattice  energy  indicates  the  correct  structure  ~80%  of  the 
time.  It  is  expected  that  the  predictive  power  of  these  methods  will  further  increase  after  the 
current  massive  re-parameterization  of  force  fields  by  the  Ammon  group  is  finished. 

Sources  of  failure  of  ab  initio  crystal  prediction  were  identified: 

1 .  Failure  to  identify  the  global  minimum  (ranking  of  candidate  crystals) 

a.  Inadequate  configuration  space  search 

h.  Rigid  molecule  approximation  when  applied  to  floppy  molecules 
c.  Lattice  energy  calculations  are  not  accurate  enough. 

2.  Kinetic  thermodynamic  problem 


Recommendations  for  addressing  these  issues  using  currently  available  computational  tools 

Recommendation  to  address  issue  l.a 

•  This  problem  has  become  less  important  as  modem  lattice  optimization  programs  can 
sample  huge  numbers  of  structures  within  reasonable  computer  time.  Although  some  unusually 
complex  crystal  stmctures  (Z'  >  1 ,  disordered,  etc.)  caimot  be  explicitly  and  accurately  found 
with  current  searches,  the  searches  are  likely  to  generate  simpler,  closely  related  stmctures  of 
comparable  energy. 

Recommendation  to  address  issues  1  .b 

•  Perform  local-density  approximation  (LDA),  since  LDA  is  known  to  give  accurate 
geometries,  or  ab  initio  calculations  of  various  conformers  of  isolated  molecule.  This  would 
include  geometry  optimizations  for  stable  stmctures  and  a  few  single-point  calculations  for 
stmctures  representing  intermediate  geometries  for  inter-conversion  from  one  conformer  to 
another.  The  number  of  conformers/stmctures  to  be  considered  depends  on  how  high  are 
conformers’  energies  above  the  global  minimum  compared  to  the  lattice  energy  of  a  given 
crystal.  LDA  is  proposed  here  as  a  very  fast  method  since  hundreds  of  calculations  during  the 
optimization  of  geometries  of  large  molecules  including  conformers  will  have  to  be  performed. 
However,  LDA  energies  are  rather  inaccurate. 

•  Use  the  molecular  stmctures  of  the  isolated  conformers  obtained  in  this  way  to  generate 
candidate  crystal  stmctures  (e.g.,  MOLPAK  program  search)  and  subsequent  energy 
minimization  assuming  simple  analytic  potential  energy  function. 
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•  Rank  candidates  in  energy  using  both  lattice  energies  and  energies  of  isolated  molecules. 

•  Perform  ab  initio  monomer  calculations  at  LDA  geometries  for  a  number  of  most 
promising  candidates  and  refine  the  energy  rankings  by  replacing  LDA  monomer  energies  by  ab 
initio  ones.  Furthermore,  ab  initio  full  geometry  optimization  can  be  performed  for  a  couple  of 
top  candidates  that  represent  local  minima  to  refine  also  monomer  geometries,  and  the  MOLPAK 
search  can  then  be  repeated. 

•  Alternatively,  periodic  boimdary  conditions  (PBC)  Moeller-Plesset  perturbation  theory  in 
second  order  (MP2)  calculations  can  be  performed  for  a  few  most  promising  candidates. 

Alternative  recommendation  to  address  issue  l.b 

•  Investigate  effect  on  packing  due  to  inclusion  of  limited  flexibility  (i.e.,  flexibility  in  only 
a  few  degrees  of  freedom  chosen  based  on  chemical  structure)  of  the  “building  block”  molecule  - 
A  flexible  force  field  as  a  function  of  small  Cartesian  displacements  of  each  atom  from  its 
equilibrium  position  can  be  obtained  from  the  energy  second  derivative  matrix  calculated  from 
isolated  molecule  calculations.  This  will  be  able  to  reflect  the  small  deviations  in  conformation, 
such  as  distortions  of  covalent  bonds  to  hydrogen  (e.g.  making  -NH2  pyramidal)  or  rotation  of 
nitro  groups,  which  are  induced  by  the  crystal  packing  to  improve  hydrogen  bonding  energies  or 
allow  denser  packings.  The  effects  of  flexibility  are  also  under  development  by  the  Ammon 
group.  Their  approach  is  to  use  intramolecular  energy  differences  (from  PM3  semi-empirical  or 
Dreiding  molecular  mechanics  calculations)  plus  intermolecular  energies  to  provide  a  total 
energy  for  the  packing  process.  One  aspect  not  being  addressed  at  present  is  the  refinement  of 
crude  structures  by  minimization  of  the  lattice  plus  intramolecular  energies.  This  will  require  at 
least  significant  modification  to  existing  software  and  may  require  development  of  a  new  suite  of 
programs. 

Recommendation  to  address  issue  1  .c 

•  Improve  the  description  of  van  der  Waals  interactions  by  using  better  asymptotic 
constants  obtained  from  ab  initio  calculations. 

•  Get  experimental  sublimation  energies  in  order  to  compare  experiment  with  calculations 
and  improve  empirical  potentials,  in  particular  for  multi-nitro  and  multi-nitramino  compounds. 

Recommendation  to  address  issue  2 

•  The  existence  of  polymorphism  shows  that  thermodynamically  metastable  crystal 
structures  are  observed.  Investigate  methods  that  the  pharmaceutical  industry  uses  to 
experimentally  identify  possible  polymorphs  and  assess  their  relative  thermodynamic  stability. 
(They  are  looking  to  computational  predictions  as  back  up  for  these  searches,  to  help  increase 
confidence  that  all  long-lived  polymorphs  have  been  found.  There  are  cases  when  the  more  stable 
polymorph  has  only  been  found  late  in  the  development  process,  and  a  few  disasters  when  it  has 
been  found  in  the  manufacturing  process). 

•  Assess  effect  of  free  energy  using  quasi-harmonic  calculation;  do  free  energy 
minimization;  if  unstable,  do  Molecular  Dynamics  on  candidate  crystals.  However,  the 
consensus  was  that  considerations  of  free  energies  make  little  difference  from  the  point  of  view 
of  crystal  structure  predictions,  as  this  only  reshuffles  the  relative  stability  when  there  are  many 
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possible  crystal  structures  within  the  energy  gap  of  possible  polymorphism.  It  does,  however, 
improve  the  realism  of  the  thermodynamic  stability  comparison. 

Recommendation  on  the  use  of  PBC  codes: 

•  Current  DFT  implementations  are  not  reliable  for  describing  the  interactions  between 
molecules  within  an  organic  crystal.  Besides,  energy  minimizations  needed  to  find  the  global 
minimum  would  be  too  time-consuming  with  these  methods. 

•  MP2  will  give  reliable  predictions  but  is  more  expensive  than  DFT.  As  stated  above,  MP2 
can  be  used  to  compute  lattice  energies  for  a  number  of  close-lying  structures  obtained  by  other 
methods.  Unfortunately,  the  size  of  basis  set  and  level  of  theory  required  to  properly  calculate  the 
dispersion  energy  will  make  this  a  considerable  challenge.  Underestimating  the  dispersion  energy 
will  tend  to  overestimate  cell  volumes,  which  will  give  errors  that  are  unlikely  to  be  comparable 
between  different  hypothetical  structures. 

Recommendations  for  futnre  research  to  advance  ab  initio  crystal  prediction 

Recommendation  to  address  issue  l.c 

•  Increase  accmacy  of  atom-atom  potentials  using  symmetry-adapted  perturbation  theory 
(SAPT)  of  intermolecular  interactions. 

•  Develop  more  realistic  potential  energy  descriptions  to  include  multi-body  effects  for  use 
in  packing  calculations. 

•  Investigate  transferability  of  the  potentials. 

•  Couple  DFT,  tight-binding  or  MNDO  code  with  the  packing  code  to  use  self-consistent 
electrostatics  in  the  solid.  This  would  allow  an  improved  electrostatic  model  that  was  consistent 
with  both  the  internal  configuration  of  the  molecule  as  well  as  the  crystalline  environment. 

•  Make  problems  available  to  people  who  are  working  on  problem  of  better  descriptions  of 
dispersion. 

•  Develop  gradients  for  MP2/PBC. 

•  Develop  DFT  functionals  that  can  be  used  for  molecules  of  interest  in  PBC  crystal  lattice 
calculations  with  plane  wave  or  atomic  basis  sets. 

Recommendation  to  address  issue  2 

•  Investigate  ways  to  simulate  nucleation/crystal  growth  to  establish  which  energetically 
feasible  crystal  structures  are  kinetically  favored.  Obtain  experimental  information  to  validate 
methods 
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Workshop  on 

Tools  for  Condensed  Phase  Computational 

Chemistry 

Clayton  Hall,  University  of  Delaware,  27-30  May  2003 

PROGRAM 


May  27:  Registration  starts  at  18:00  and  dinner  at  19:00  in  Clayton  Haii. 

May  28:  Registration  7:30-10:00  in  Ciayton  Hall. 

DOD  Programs  and  Needs  for  Computational  Chemistry  Packages 


08:30-09:15 

Betsy  Rice,  Cary  Chaba-  - 
lowski,  William  Mattson 

ARL 

Army  Needs  for  Crystal  Structure  and  other  Property 
Predictions 

09:15-09:45 

Betsy  Rice 

ARL 

Parallelizing  Molpak  and  Planewave 

09:45-10:00 

Break 

10:00-10:30 

David  Singh 

NRL 

Pseudopotential  and  Planewave  Codes 

Presentations  by  Experts 


10:30-11:30 

Krzysztof  Szalewicz 

U.  Delaware 

Intermolecular  Interaction  Potentials 

11:30-12:30 

Herman  Ammon 

U.  Maryland 

Crystal  Structure  Prediction  for  Energetic  Materials 

12:30-13:30 

Lunch 

13:30:14:30 

Sally  Price 

U.  Coll. 

London 

Anisotropic  Atom-Atom  Intermolecular  Potentials  in 
Organic  Crystal  Structure  and  Property  Prediction 

14:30:15:30 

Emilio  Artacho 

Cambridge  U. 

Material  Properties  with  SIESTA:  Strengths, 
Weaknesses,  and  Prospects 

15:30-16:00 

Break 

16::00-17:00 

Julian  Gale 

Imperial  Coll. 

General  Utility  Lattice  Program 

17:00-18:00 

Anne  Chaka 

NIST 

Validation  of  Computational  Methods 

18:00-19:30 

Dinner 

19:30-21:00 

Discussion 

Possible  Implementations  of  New  Approaches  by 

Army 
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May  29: 

Presentations  by  Experts  (continued) 


0830-09:30 

Eric  Bylaska 

PNNL 

Materials  Properties  with  NWChem 

09:30-10:30 

Richard  Martin 

U.  Illinois 

Identifying  the  Key  Problems  in  Present  Density 
Functionals 

10:30-11:00 

Break 

11:00-12:00 

Gustavo  Scuseria 

Rice  U. 

Condensed  Phase  Simulations  using  Gaussian 

Orbitals  and  Periodic  Boundary  Conditions 

12:00-13:00 

Lunch 

13:00-14:00 

Emily  Carter 

UCLA 

Modeling  Chemistry  and  Physics  in  Bulk  Crystals, 
Surfaces,  and  Interfaces:  What's  Possible  and  What 
Can  We  Trust? 

Group  discussion:  Best  Ways  to  Proceed  on  Army  Probiems 


14:00 

Krzysztof  Szalewicz  and 
Bob  Shaw 

Summary  remarks  and  guidelines  for  discussion 
groups 

14:00-17:00 

Discussion  group  1 

Chabalowski,  Byrd,  Mattson,  Martin,  Carter,  Bylaska, 
Scuseria,  Doren,  Tchoukova,  Murdachaew,  Sandler, 
Shaw 

14:00-17:00 

Discussion  group  II 

Rice,  Kim,  Ammon,  Price,  Gale,  Artacho,  Chaka, 
Singh,  Bukowski,  Akin-Ojo,  Wang,  Rowe,  Szalewicz 

15:30-15:45 

Break 

17:00-18:00 

Discussion 

Groups  Report  back  with  their  Conclusions 

18:00 

Dinner 
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